# Replication code for "What is rebel governance? Introducing a new 
# dataset on rebel institutions, 1945–2012"
# 
# Original author: Karen Albert
# Replicated by: Daniel Korell, Pelke, Lars
# R version R-4.2.0
#----------------------------------------------------------------------------- #

library(broom) # Data tidying
library(tidyverse) # Data manipulation & create figures
library(modelsummary) # Create figures and tables
library(lfe) # Fixed Effects
library(fixest)
library(patchwork)

# Set a working directory to where the data is saved
setwd("C:/Users/ba72loko/projects/replication paper ZfVP/calculations") 

rm(list = ls())

# Load data via csv
qsi <- read.csv('data/Study_1/RQSI_Replication_data_2021.csv', stringsAsFactors = FALSE)

# load V-Dem data (v12) #
vdem.data <- readRDS("data/Study_1/V-Dem-CY-Full+Others-v12.rds")

# merge datasets#

qsi.vdem <- qsi %>%
  dplyr::select(-X.1) %>%
  rename(country_id = VDemID, 
         year = Year) %>%
  left_join(vdem.data, by = c("country_id", "year")) %>%
  mutate(e_pop_log = log10(e_pop+1), 
         e_total_resources_income_pc_log = log10(e_total_resources_income_pc+1))

################################################################################

# Replication of the original study estimators

################################################################################

# Fixed Effects - Year
# Rebel institutions and relative rebel strength

dv <- c("Constitution.count.1", "Law.count.1", "JoinIO.count.1", "DipMiss.count.1", "BorderP.count.1", "Justice.count.1", 
        "Police.count.1", "Elections.count.1", "Currency.count.1", "ParallelGov.count.1", "LocalGov.count.1", 
        "Gov_in_Exile.count.1", "Tax.count.1", "EconTreaty.count.1", "NegResource.count.1", "Edu.count.1", "Health.count.1")

for(i in 1:length(dv)){
  
  model <- paste("m",i, sep="_")
  
  m <- feols(as.formula(paste(dv[i], "~ rebstrength.num + grptime2")), data = qsi.vdem)
  
  assign(model,m)}

# Rebel institutions and technology of rebellion

for(i in 1:length(dv)){
  
  model <- paste("m_a",i, sep="_")
  
  m <- feols(as.formula(paste(dv[i], "~ technologyrebellion + grptime2")), data = qsi.vdem)
  
  assign(model,m)}



modelsummary::modelsummary(list(m_1, m_2, m_3, m_4, m_5, m_6, m_7, m_8, m_9, m_10, 
                                m_11, m_12, m_13, m_14, m_15, m_16, m_17), 
                           output = "tables/Study_1/Table_A1_original_estimates.docx", 
                           estimate = "{estimate}{stars}", 
                           statistic = 'std.error',
                           coef_rename = c("rebstrength.num" = "Rebel Strength", 
                                           "grptime2" = "Group existance count"))



# Figures

# Create a background reference line for all figures
b <- list(geom_vline(xintercept = 0, color = 'red'))

# Visualizing effect sizes

# Rebel strength

f1_data_orig <-  modelplot(list("Constitution" = m_1,
                           "Law" = m_2,
                           "Join IO" = m_3,
                           "Diplomatic mission" = m_4,
                           "Border patrol" = m_5,
                           "Justice" = m_6,
                           "Policing" = m_7,
                           "Elections" = m_8,
                           "Currency" = m_9,
                           "Parallel government" = m_10,
                           "Local government" = m_11,
                           "Government in exile" = m_12,
                           "Tax" = m_13,
                           "Economic treaty" = m_14,
                           "Negotiate resource rights" = m_15,
                           "Education" = m_16,
                           "Healthcare" = m_17),
                      coef_rename = c("technologyrebellionIrregular" = "Irregular", 
                                      "technologyrebellionSNC" = "SNC", 
                                      "grptime2" = "Group existance count", 
                                      "rebstrength.num" = "Relative rebel strength",
                                      "e_peinfmor" = "Infant mortality",
                                      "v2exl_legitideol" = "Ideology",
                                      "e_total_resources_income_pc_log" = "Resource income",
                                      "e_pop_log" = "Population",
                                      "e_gdppc_log" = "GDP per capita"),
                      coef_omit = 'Intercept|e_peinfmor|v2exl_legitideol|
          |e_total_resources_income_pc_log|e_gdppc|technologyrebellionSNC|
          |e_pop_log|grptime2', 
                      facet=T, background = b,
                      draw = FALSE)


# Technology of the Rebellion

f2_data_orig <- modelplot(list("Constitution" = m_a_1,
                          "Law" = m_a_2,
                          "Join IO" = m_a_4,
                          "Diplomatic mission" = m_a_4,
                          "Border patrol" = m_a_5,
                          "Justice" = m_a_6,
                          "Policing" = m_a_7,
                          "Elections" = m_a_8,
                          "Currency" = m_a_9,
                          "Parallel government" = m_a_10,
                          "Local government" = m_a_11,
                          "Government in exile" = m_a_12,
                          "Tax" = m_a_13,
                          "Economic treaty" = m_a_14,
                          "Negotiate resource rights" = m_a_15,
                          "Education" = m_a_16,
                          "Healthcare" = m_a_17), 
                     coef_rename = c("technologyrebellionIrregular" = "Irregular", 
                                     "technologyrebellionSNC" = "SNC", 
                                     "grptime2" = "Group existance count", 
                                     "rebstrength.num" = "Relative rebel strength",
                                     "e_peinfmor" = "Infant mortality",
                                     "v2exl_legitideol" = "Ideology",
                                     "e_total_resources_income_pc_log" = "Resource income",
                                     "e_pop_log" = "Population",
                                     "e_gdppc_log" = "GDP per capita"),
                     coef_omit = 'Intercept|e_peinfmor|v2exl_legitideol|
          |e_total_resources_income_pc_log|e_gdppc|technologyrebellionSNC|
          |e_pop_log|grptime2', 
                     facet=T, background = b, draw = F) 

theme_set(theme_minimal())

# plot
f1_orig <- ggplot(f1_data_orig, aes(
  y = model, x = estimate,
  xmin = conf.low, xmax = conf.high)) +
  ggdist::geom_pointinterval() +
  geom_vline(xintercept = 0, color = 'red', linetype = "longdash",  alpha = 0.5) +
  xlab("Estimate\nRelative Rebel Strength") +
  ylab("") +
  ggtitle("Original estimates")


# plot
f2_orig <- ggplot(f2_data_orig, aes(
  y = model, x = estimate,
  xmin = conf.low, xmax = conf.high)) +
  ggdist::geom_pointinterval() +
  geom_vline(xintercept = 0, color = 'red', linetype = "longdash", alpha = 0.5) +
  xlab("Estimate\nTechnology of Rebellion") +
  ylab("") +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank()) 


################################################################################

# Main Analysis: Regression Models Main Paper with Two-Way-Fixed Effects and Control Variables

################################################################################

# Fixed Effects - Year
# Rebel institutions and relative rebel strength

dv <- c("Constitution.count.1", "Law.count.1", "JoinIO.count.1", "DipMiss.count.1", "BorderP.count.1", "Justice.count.1", 
        "Police.count.1", "Elections.count.1", "Currency.count.1", "ParallelGov.count.1", "LocalGov.count.1", 
        "Gov_in_Exile.count.1", "Tax.count.1", "EconTreaty.count.1", "NegResource.count.1", "Edu.count.1", "Health.count.1")

for(i in 1:length(dv)){
  
  model <- paste("m",i, "twfe", sep="_")
  
  m <- feols(as.formula(paste(dv[i], "~ rebstrength.num + e_gdppc + 
                  e_pop_log + e_total_resources_income_pc_log + v2exl_legitideol + 
                              e_peinfmor | Location + year")), cluster = "Location",
             data = qsi.vdem)
  
  assign(model,m)}

# Rebel institutions and technology of rebellion

for(i in 1:length(dv)){
  
  model <- paste("m_a",i, "twfe", sep="_")
  
  m <- feols(as.formula(paste(dv[i], "~ technologyrebellion + e_gdppc + 
                  e_pop_log + e_total_resources_income_pc_log + v2exl_legitideol + 
                              e_peinfmor | Location + year")), cluster = "Location",
             data = qsi.vdem)
  
  assign(model,m)}

modelsummary::modelsummary(list(m_1_twfe, m_2_twfe, m_3_twfe, m_4_twfe, m_5_twfe, m_6_twfe, m_7_twfe, m_8_twfe, m_9_twfe, m_10_twfe, 
                                m_11_twfe, m_12_twfe, m_13_twfe, m_14_twfe, m_15_twfe, m_16_twfe, m_17_twfe), 
                           output = "tables/Study_1/Table_A1.docx", 
                           estimate = "{estimate}{stars}", 
                           statistic = 'std.error',
                           coef_rename = c("rebstrength.num" = "Rebel Strength", 
                                           "e_gdppc" = "GDP pc (log10)", 
                                           "e_pop_log" = "Pop (log10)", 
                                           "e_total_resources_income_pc_log" = "Ressource Income (log10)", 
                                           "v2exl_legitideol" = "Ideology", 
                                           "e_peinfmor" = "Infant Mortality"))



# Figures

# Create a background reference line for all figures
b <- list(geom_vline(xintercept = 0, color = 'red'))

# Visualizing effect sizes

# Rebel strength

f1_data <-  modelplot(list("Constitution" = m_1_twfe,
                           "Law" = m_2_twfe,
                           "Join IO" = m_3_twfe,
                           "Diplomatic mission" = m_4_twfe,
                           "Border patrol" = m_5_twfe,
                           "Justice" = m_6_twfe,
                           "Policing" = m_7_twfe,
                           "Elections" = m_8_twfe,
                           "Currency" = m_9_twfe,
                           "Parallel government" = m_10_twfe,
                           "Local government" = m_11_twfe,
                           "Government in exile" = m_12_twfe,
                           "Tax" = m_13_twfe,
                           "Economic treaty" = m_14_twfe,
                           "Negotiate resource rights" = m_15_twfe,
                           "Education" = m_16_twfe,
                           "Healthcare" = m_17_twfe),
                      coef_rename = c("technologyrebellionIrregular" = "Irregular", 
                                      "technologyrebellionSNC" = "SNC", 
                                      "grptime2" = "Group existance count", 
                                      "rebstrength.num" = "Relative rebel strength",
                                      "e_peinfmor" = "Infant mortality",
                                      "v2exl_legitideol" = "Ideology",
                                      "e_total_resources_income_pc_log" = "Resource income",
                                      "e_pop_log" = "Population",
                                      "e_gdppc_log" = "GDP per capita"),
                      coef_omit = 'Intercept|e_peinfmor|v2exl_legitideol|
          |e_total_resources_income_pc_log|e_gdppc|technologyrebellionSNC|
          |e_pop_log', 
                      facet=T, background = b,
                      draw = FALSE)


# Technology of the Rebellion

f2_data <- modelplot(list("Constitution" = m_a_1_twfe,
               "Law" = m_a_2_twfe,
               "Join IO" = m_a_4_twfe,
               "Diplomatic mission" = m_a_4_twfe,
               "Border patrol" = m_a_5_twfe,
               "Justice" = m_a_6_twfe,
               "Policing" = m_a_7_twfe,
               "Elections" = m_a_8_twfe,
               "Currency" = m_a_9_twfe,
               "Parallel government" = m_a_10_twfe,
               "Local government" = m_a_11_twfe,
               "Government in exile" = m_a_12_twfe,
               "Tax" = m_a_13_twfe,
               "Economic treaty" = m_a_14_twfe,
               "Negotiate resource rights" = m_a_15_twfe,
               "Education" = m_a_16_twfe,
               "Healthcare" = m_a_17_twfe), 
          coef_rename = c("technologyrebellionIrregular" = "Irregular", 
                          "technologyrebellionSNC" = "SNC", 
                          "grptime2" = "Group existance count", 
                          "rebstrength.num" = "Relative rebel strength",
                          "e_peinfmor" = "Infant mortality",
                          "v2exl_legitideol" = "Ideology",
                          "e_total_resources_income_pc_log" = "Resource income",
                          "e_pop_log" = "Population",
                          "e_gdppc_log" = "GDP per capita"),
          coef_omit = 'Intercept|e_peinfmor|v2exl_legitideol|
          |e_total_resources_income_pc_log|e_gdppc|technologyrebellionSNC|
          |e_pop_log', 
          facet=T, background = b, draw = F) 

theme_set(theme_minimal())

# plot
f1 <- ggplot(f1_data, aes(
  y = model, x = estimate,
  xmin = conf.low, xmax = conf.high)) +
  ggdist::geom_pointinterval() +
  geom_vline(xintercept = 0, color = 'red', linetype = "longdash",  alpha = 0.5) +
  xlab("Estimate\nRelative Rebel Strength") +
  ylab("") +
  ggtitle("Replication with TWFE and controls")
  
# plot
f2 <- ggplot(f2_data, aes(
  y = model, x = estimate,
  xmin = conf.low, xmax = conf.high)) +
  ggdist::geom_pointinterval() +
  geom_vline(xintercept = 0, color = 'red', linetype = "longdash", alpha = 0.5) +
    xlab("Estimate\nTechnology of Rebellion") +
  ylab("") +
  theme(axis.title.y = element_blank(),
                 axis.text.y = element_blank(),
                 axis.ticks.y = element_blank()) 

(f1_orig + f2_orig)/ (f1 + f2) +  plot_annotation(tag_levels = 'A')

ggsave("figures/Study_1/Figure_1_Main_paper.png", width = 17, height = 20, units = "cm", dpi = 600)



################################################################################

# Appendix: Regression Models with Region Fixed Effects

################################################################################

# Fixed Effects - Year
# Rebel institutions and relative rebel strength

dv <- c("Constitution.count.1", "Law.count.1", "JoinIO.count.1", "DipMiss.count.1", "BorderP.count.1", "Justice.count.1", 
        "Police.count.1", "Elections.count.1", "Currency.count.1", "ParallelGov.count.1", "LocalGov.count.1", 
        "Gov_in_Exile.count.1", "Tax.count.1", "EconTreaty.count.1", "NegResource.count.1", "Edu.count.1", "Health.count.1")

for(i in 1:length(dv)){
  
  model <- paste("m",i, "twfe", sep="_")
  
  m <- feols(as.formula(paste(dv[i], "~ rebstrength.num + e_gdppc + 
                  e_pop_log + e_total_resources_income_pc_log + v2exl_legitideol + 
                              e_peinfmor | region + year")), cluster = "Location",
             data = qsi.vdem)
  
  assign(model,m)}

# Rebel institutions and technology of rebellion

for(i in 1:length(dv)){
  
  model <- paste("m_a",i, "twfe", sep="_")
  
  m <- feols(as.formula(paste(dv[i], "~ technologyrebellion + e_gdppc + 
                  e_pop_log + e_total_resources_income_pc_log + v2exl_legitideol + 
                              e_peinfmor | region + year")), cluster = "Location",
             data = qsi.vdem)
  
  assign(model,m)}

modelsummary::modelsummary(list(m_1_twfe, m_2_twfe, m_3_twfe, m_4_twfe, m_5_twfe, m_6_twfe, m_7_twfe, m_8_twfe, m_9_twfe, m_10_twfe, 
                                m_11_twfe, m_12_twfe, m_13_twfe, m_14_twfe, m_15_twfe, m_16_twfe, m_17_twfe), 
                           output = "tables/Study_1/Table_A2_Region_FE.docx", 
                           estimate = "{estimate}{stars}", 
                           statistic = 'std.error',
                           coef_rename = c("rebstrength.num" = "Rebel Strength", 
                                           "e_gdppc" = "GDP pc (log10)", 
                                           "e_pop_log" = "Pop (log10)", 
                                           "e_total_resources_income_pc_log" = "Ressource Income (log10)", 
                                           "v2exl_legitideol" = "Ideology", 
                                           "e_peinfmor" = "Infant Mortality"))



# Figures

# Create a background reference line for all figures
b <- list(geom_vline(xintercept = 0, color = 'red'))

# Visualizing effect sizes

# Rebel strength

f1_data <-  modelplot(list("Constitution" = m_1_twfe,
                           "Law" = m_2_twfe,
                           "Join IO" = m_3_twfe,
                           "Diplomatic mission" = m_4_twfe,
                           "Border patrol" = m_5_twfe,
                           "Justice" = m_6_twfe,
                           "Policing" = m_7_twfe,
                           "Elections" = m_8_twfe,
                           "Currency" = m_9_twfe,
                           "Parallel government" = m_10_twfe,
                           "Local government" = m_11_twfe,
                           "Government in exile" = m_12_twfe,
                           "Tax" = m_13_twfe,
                           "Economic treaty" = m_14_twfe,
                           "Negotiate resource rights" = m_15_twfe,
                           "Education" = m_16_twfe,
                           "Healthcare" = m_17_twfe),
                      coef_rename = c("technologyrebellionIrregular" = "Irregular", 
                                      "technologyrebellionSNC" = "SNC", 
                                      "grptime2" = "Group existance count", 
                                      "rebstrength.num" = "Relative rebel strength",
                                      "e_peinfmor" = "Infant mortality",
                                      "v2exl_legitideol" = "Ideology",
                                      "e_total_resources_income_pc_log" = "Resource income",
                                      "e_pop_log" = "Population",
                                      "e_gdppc_log" = "GDP per capita"),
                      coef_omit = 'Intercept|e_peinfmor|v2exl_legitideol|
          |e_total_resources_income_pc_log|e_gdppc|technologyrebellionSNC|
          |e_pop_log', 
                      facet=T, background = b,
                      draw = FALSE)


# Technology of the Rebellion

f2_data <- modelplot(list("Constitution" = m_a_1_twfe,
                          "Law" = m_a_2_twfe,
                          "Join IO" = m_a_4_twfe,
                          "Diplomatic mission" = m_a_4_twfe,
                          "Border patrol" = m_a_5_twfe,
                          "Justice" = m_a_6_twfe,
                          "Policing" = m_a_7_twfe,
                          "Elections" = m_a_8_twfe,
                          "Currency" = m_a_9_twfe,
                          "Parallel government" = m_a_10_twfe,
                          "Local government" = m_a_11_twfe,
                          "Government in exile" = m_a_12_twfe,
                          "Tax" = m_a_13_twfe,
                          "Economic treaty" = m_a_14_twfe,
                          "Negotiate resource rights" = m_a_15_twfe,
                          "Education" = m_a_16_twfe,
                          "Healthcare" = m_a_17_twfe), 
                     coef_rename = c("technologyrebellionIrregular" = "Irregular", 
                                     "technologyrebellionSNC" = "SNC", 
                                     "grptime2" = "Group existance count", 
                                     "rebstrength.num" = "Relative rebel strength",
                                     "e_peinfmor" = "Infant mortality",
                                     "v2exl_legitideol" = "Ideology",
                                     "e_total_resources_income_pc_log" = "Resource income",
                                     "e_pop_log" = "Population",
                                     "e_gdppc_log" = "GDP per capita"),
                     coef_omit = 'Intercept|e_peinfmor|v2exl_legitideol|
          |e_total_resources_income_pc_log|e_gdppc|technologyrebellionSNC|
          |e_pop_log', 
                     facet=T, background = b, draw = F) 


# plot
f1 <- ggplot(f1_data, aes(
  y = model, x = estimate,
  xmin = conf.low, xmax = conf.high)) +
  ggdist::geom_pointinterval() +
  geom_vline(xintercept = 0, color = 'red', linetype = "longdash",  alpha = 0.5) +
  xlab("Estimate") +
  ylab("Dependent Variable") 

# plot
f2 <- ggplot(f2_data, aes(
  y = model, x = estimate,
  xmin = conf.low, xmax = conf.high)) +
  ggdist::geom_pointinterval() +
  geom_vline(xintercept = 0, color = 'red', linetype = "longdash", alpha = 0.5) +
  xlab("Estimate") +
  ylab("") +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank())

f1 + f2
ggsave("figures/Study_1/Figure_2_Region_FE.png", width = 15, height = 10, units = "cm", dpi = 600)


################################################################################

# Appendix: Regression Models for Asia

################################################################################

# Fixed Effects - Year
# Rebel institutions and relative rebel strength

dv <- c("Constitution.count.1", "Law.count.1", "JoinIO.count.1", "DipMiss.count.1", "BorderP.count.1", "Justice.count.1", 
        "Police.count.1", "Elections.count.1", "Currency.count.1", "ParallelGov.count.1", "LocalGov.count.1", 
        "Gov_in_Exile.count.1", "Tax.count.1", "EconTreaty.count.1", "NegResource.count.1", "Edu.count.1", "Health.count.1")

for(i in 1:length(dv)){
  
  model <- paste("m",i, "twfe", sep="_")
  
  m <- feols(as.formula(paste(dv[i], "~ rebstrength.num + e_gdppc + 
                  e_pop_log + e_total_resources_income_pc_log + v2exl_legitideol + 
                              e_peinfmor")), cluster = "Location",
             data = subset(qsi.vdem, region=="Asia"))
  
  assign(model,m)}

# Rebel institutions and technology of rebellion

for(i in 1:length(dv)){
  
  model <- paste("m_a",i, "twfe", sep="_")
  
  m <- feols(as.formula(paste(dv[i], "~ technologyrebellion + e_gdppc + 
                  e_pop_log + e_total_resources_income_pc_log + v2exl_legitideol + 
                              e_peinfmor")), cluster = "Location",
             data = subset(qsi.vdem, region=="Asia"))
  
  assign(model,m)}

modelsummary::modelsummary(list(m_1_twfe, m_2_twfe, m_3_twfe, m_4_twfe, m_5_twfe, m_6_twfe, m_7_twfe, m_8_twfe, m_9_twfe, m_10_twfe, 
                                m_11_twfe, m_12_twfe, m_13_twfe, m_14_twfe, m_15_twfe, m_16_twfe, m_17_twfe), 
                           output = "tables/Study_1/Table_A3_Asia.docx", 
                           estimate = "{estimate}{stars}", 
                           statistic = 'std.error',
                           coef_rename = c("rebstrength.num" = "Rebel Strength", 
                                           "e_gdppc" = "GDP pc (log10)", 
                                           "e_pop_log" = "Pop (log10)", 
                                           "e_total_resources_income_pc_log" = "Ressource Income (log10)", 
                                           "v2exl_legitideol" = "Ideology", 
                                           "e_peinfmor" = "Infant Mortality"))



# Figures

# Create a background reference line for all figures
b <- list(geom_vline(xintercept = 0, color = 'red'))

# Visualizing effect sizes

# Rebel strength

f1_data <-  modelplot(list("Constitution" = m_1_twfe,
                           "Law" = m_2_twfe,
                           "Join IO" = m_3_twfe,
                           "Diplomatic mission" = m_4_twfe,
                           "Border patrol" = m_5_twfe,
                           "Justice" = m_6_twfe,
                           "Policing" = m_7_twfe,
                           "Elections" = m_8_twfe,
                           "Currency" = m_9_twfe,
                           "Parallel government" = m_10_twfe,
                           "Local government" = m_11_twfe,
                           "Government in exile" = m_12_twfe,
                           "Tax" = m_13_twfe,
                           "Economic treaty" = m_14_twfe,
                           "Negotiate resource rights" = m_15_twfe,
                           "Education" = m_16_twfe,
                           "Healthcare" = m_17_twfe),
                      coef_rename = c("technologyrebellionIrregular" = "Irregular", 
                                      "technologyrebellionSNC" = "SNC", 
                                      "grptime2" = "Group existance count", 
                                      "rebstrength.num" = "Relative rebel strength",
                                      "e_peinfmor" = "Infant mortality",
                                      "v2exl_legitideol" = "Ideology",
                                      "e_total_resources_income_pc_log" = "Resource income",
                                      "e_pop_log" = "Population",
                                      "e_gdppc_log" = "GDP per capita"),
                      coef_omit = 'Intercept|e_peinfmor|v2exl_legitideol|
          |e_total_resources_income_pc_log|e_gdppc|technologyrebellionSNC|
          |e_pop_log', 
                      facet=T, background = b,
                      draw = FALSE)


# Technology of the Rebellion

f2_data <- modelplot(list("Constitution" = m_a_1_twfe,
                          "Law" = m_a_2_twfe,
                          "Join IO" = m_a_4_twfe,
                          "Diplomatic mission" = m_a_4_twfe,
                          "Border patrol" = m_a_5_twfe,
                          "Justice" = m_a_6_twfe,
                          "Policing" = m_a_7_twfe,
                          "Elections" = m_a_8_twfe,
                          "Currency" = m_a_9_twfe,
                          "Parallel government" = m_a_10_twfe,
                          "Local government" = m_a_11_twfe,
                          "Government in exile" = m_a_12_twfe,
                          "Tax" = m_a_13_twfe,
                          "Economic treaty" = m_a_14_twfe,
                          "Negotiate resource rights" = m_a_15_twfe,
                          "Education" = m_a_16_twfe,
                          "Healthcare" = m_a_17_twfe), 
                     coef_rename = c("technologyrebellionIrregular" = "Irregular", 
                                     "technologyrebellionSNC" = "SNC", 
                                     "grptime2" = "Group existance count", 
                                     "rebstrength.num" = "Relative rebel strength",
                                     "e_peinfmor" = "Infant mortality",
                                     "v2exl_legitideol" = "Ideology",
                                     "e_total_resources_income_pc_log" = "Resource income",
                                     "e_pop_log" = "Population",
                                     "e_gdppc_log" = "GDP per capita"),
                     coef_omit = 'Intercept|e_peinfmor|v2exl_legitideol|
          |e_total_resources_income_pc_log|e_gdppc|technologyrebellionSNC|
          |e_pop_log', 
                     facet=T, background = b, draw = F) 


# plot
f1 <- ggplot(f1_data, aes(
  y = model, x = estimate,
  xmin = conf.low, xmax = conf.high)) +
  ggdist::geom_pointinterval() +
  geom_vline(xintercept = 0, color = 'red', linetype = "longdash",  alpha = 0.5) +
  xlab("Estimate") +
  ylab("Dependent Variable") 

# plot
f2 <- ggplot(f2_data, aes(
  y = model, x = estimate,
  xmin = conf.low, xmax = conf.high)) +
  ggdist::geom_pointinterval() +
  geom_vline(xintercept = 0, color = 'red', linetype = "longdash", alpha = 0.5) +
  xlab("Estimate") +
  ylab("") +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank())

f1 + f2
ggsave("figures/Study_1/Figure_3_Asia.png", width = 15, height = 10, units = "cm", dpi = 600)



################################################################################

# Appendix: Regression Models for Sub-Sahara-Africa

################################################################################

# Fixed Effects - Year
# Rebel institutions and relative rebel strength

dv <- c("Constitution.count.1", "Law.count.1", "JoinIO.count.1", "DipMiss.count.1", "BorderP.count.1", "Justice.count.1", 
        "Police.count.1", "Elections.count.1", "Currency.count.1", "ParallelGov.count.1", "LocalGov.count.1", 
        "Gov_in_Exile.count.1", "Tax.count.1", "EconTreaty.count.1", "NegResource.count.1", "Edu.count.1", "Health.count.1")

for(i in 1:length(dv)){
  
  model <- paste("m",i, "twfe", sep="_")
  
  m <- feols(as.formula(paste(dv[i], "~ rebstrength.num + e_gdppc + 
                  e_pop_log + e_total_resources_income_pc_log + v2exl_legitideol + 
                              e_peinfmor | year")), cluster = "Location",
             data = subset(qsi.vdem, region=="SubAfrica"))
  
  assign(model,m)}

# Rebel institutions and technology of rebellion

for(i in 1:length(dv)){
  
  model <- paste("m_a",i, "twfe", sep="_")
  
  m <- feols(as.formula(paste(dv[i], "~ technologyrebellion + e_gdppc + 
                  e_pop_log + e_total_resources_income_pc_log + v2exl_legitideol + 
                              e_peinfmor | year")), cluster = "Location",
             data = subset(qsi.vdem, region=="SubAfrica"))
  
  assign(model,m)}

modelsummary::modelsummary(list(m_1_twfe, m_2_twfe, m_3_twfe, m_4_twfe, m_5_twfe, m_6_twfe, m_7_twfe, m_8_twfe, m_9_twfe, m_10_twfe, 
                                m_11_twfe, m_12_twfe, m_13_twfe, m_14_twfe, m_15_twfe, m_16_twfe, m_17_twfe), 
                           output = "tables/Study_1/Table_A3_SSA.docx", 
                           estimate = "{estimate}{stars}", 
                           statistic = 'std.error',
                           coef_rename = c("rebstrength.num" = "Rebel Strength", 
                                           "e_gdppc" = "GDP pc (log10)", 
                                           "e_pop_log" = "Pop (log10)", 
                                           "e_total_resources_income_pc_log" = "Ressource Income (log10)", 
                                           "v2exl_legitideol" = "Ideology", 
                                           "e_peinfmor" = "Infant Mortality"))



# Figures

# Rebel strength

f1_data <-  modelplot(list("Constitution" = m_1_twfe,
                           "Law" = m_2_twfe,
                           "Join IO" = m_3_twfe,
                           "Diplomatic mission" = m_4_twfe,
                           "Border patrol" = m_5_twfe,
                           "Justice" = m_6_twfe,
                           "Policing" = m_7_twfe,
                           "Elections" = m_8_twfe,
                           "Currency" = m_9_twfe,
                           "Parallel government" = m_10_twfe,
                           "Local government" = m_11_twfe,
                           "Government in exile" = m_12_twfe,
                           "Tax" = m_13_twfe,
                           "Economic treaty" = m_14_twfe,
                           "Negotiate resource rights" = m_15_twfe,
                           "Education" = m_16_twfe,
                           "Healthcare" = m_17_twfe),
                      coef_rename = c("technologyrebellionIrregular" = "Irregular", 
                                      "technologyrebellionSNC" = "SNC", 
                                      "grptime2" = "Group existance count", 
                                      "rebstrength.num" = "Relative rebel strength",
                                      "e_peinfmor" = "Infant mortality",
                                      "v2exl_legitideol" = "Ideology",
                                      "e_total_resources_income_pc_log" = "Resource income",
                                      "e_pop_log" = "Population",
                                      "e_gdppc_log" = "GDP per capita"),
                      coef_omit = 'Intercept|e_peinfmor|v2exl_legitideol|
          |e_total_resources_income_pc_log|e_gdppc|technologyrebellionSNC|
          |e_pop_log', 
                      facet=T, background = b,
                      draw = FALSE)


# Technology of the Rebellion

f2_data <- modelplot(list("Constitution" = m_a_1_twfe,
                          "Law" = m_a_2_twfe,
                          "Join IO" = m_a_4_twfe,
                          "Diplomatic mission" = m_a_4_twfe,
                          "Border patrol" = m_a_5_twfe,
                          "Justice" = m_a_6_twfe,
                          "Policing" = m_a_7_twfe,
                          "Elections" = m_a_8_twfe,
                          "Currency" = m_a_9_twfe,
                          "Parallel government" = m_a_10_twfe,
                          "Local government" = m_a_11_twfe,
                          "Government in exile" = m_a_12_twfe,
                          "Tax" = m_a_13_twfe,
                          "Economic treaty" = m_a_14_twfe,
                          "Negotiate resource rights" = m_a_15_twfe,
                          "Education" = m_a_16_twfe,
                          "Healthcare" = m_a_17_twfe), 
                     coef_rename = c("technologyrebellionIrregular" = "Irregular", 
                                     "technologyrebellionSNC" = "SNC", 
                                     "grptime2" = "Group existance count", 
                                     "rebstrength.num" = "Relative rebel strength",
                                     "e_peinfmor" = "Infant mortality",
                                     "v2exl_legitideol" = "Ideology",
                                     "e_total_resources_income_pc_log" = "Resource income",
                                     "e_pop_log" = "Population",
                                     "e_gdppc_log" = "GDP per capita"),
                     coef_omit = 'Intercept|e_peinfmor|v2exl_legitideol|
          |e_total_resources_income_pc_log|e_gdppc|technologyrebellionSNC|
          |e_pop_log', 
                     facet=T, background = b, draw = F) 


# plot
f1 <- ggplot(f1_data, aes(
  y = model, x = estimate,
  xmin = conf.low, xmax = conf.high)) +
  ggdist::geom_pointinterval() +
  geom_vline(xintercept = 0, color = 'red', linetype = "longdash",  alpha = 0.5) +
  xlab("Estimate") +
  ylab("Dependent Variable") 

# plot
f2 <- ggplot(f2_data, aes(
  y = model, x = estimate,
  xmin = conf.low, xmax = conf.high)) +
  ggdist::geom_pointinterval() +
  geom_vline(xintercept = 0, color = 'red', linetype = "longdash", alpha = 0.5) +
  xlab("Estimate") +
  ylab("") +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank())

f1 + f2
ggsave("figures/Study_1/Figure_3_SSA.png", width = 15, height = 10, units = "cm", dpi = 600)



################################################################################

# Appendix: Regression Models for Cold War

################################################################################

# Fixed Effects - Year
# Rebel institutions and relative rebel strength

dv <- c("Constitution.count.1", "Law.count.1", "JoinIO.count.1", "DipMiss.count.1", "BorderP.count.1", "Justice.count.1", 
        "Police.count.1", "Elections.count.1", "Currency.count.1", "ParallelGov.count.1", "LocalGov.count.1", 
        "Gov_in_Exile.count.1", "Tax.count.1", "EconTreaty.count.1", "NegResource.count.1", "Edu.count.1", "Health.count.1")

for(i in 1:length(dv)){
  
  model <- paste("m",i, "twfe", sep="_")
  
  m <- feols(as.formula(paste(dv[i], "~ rebstrength.num + e_gdppc + 
                  e_pop_log + e_total_resources_income_pc_log + v2exl_legitideol + 
                              e_peinfmor | Location + year")), cluster = "Location",
             data = subset(qsi.vdem, year %in% c(1945:1991)))
  
  assign(model,m)}

# Rebel institutions and technology of rebellion

for(i in 1:length(dv)){
  
  model <- paste("m_a",i, "twfe", sep="_")
  
  m <- feols(as.formula(paste(dv[i], "~ technologyrebellion + e_gdppc + 
                  e_pop_log + e_total_resources_income_pc_log + v2exl_legitideol + 
                              e_peinfmor | Location + year")), cluster = "Location",
             data = subset(qsi.vdem, year %in% c(1945:1991)))
  
  assign(model,m)}

modelsummary::modelsummary(list(m_1_twfe, m_2_twfe, m_3_twfe, m_4_twfe, m_5_twfe, m_6_twfe, m_7_twfe, m_8_twfe, m_9_twfe, m_10_twfe, 
                                m_11_twfe, m_12_twfe, m_13_twfe, m_14_twfe, m_15_twfe, m_16_twfe, m_17_twfe), 
                           output = "tables/Study_1/Table_A4_ColdWar.docx", 
                           estimate = "{estimate}{stars}", 
                           statistic = 'std.error',
                           coef_rename = c("rebstrength.num" = "Rebel Strength", 
                                           "e_gdppc" = "GDP pc (log10)", 
                                           "e_pop_log" = "Pop (log10)", 
                                           "e_total_resources_income_pc_log" = "Ressource Income (log10)", 
                                           "v2exl_legitideol" = "Ideology", 
                                           "e_peinfmor" = "Infant Mortality"))



# Figures

# Rebel strength

f1_data <-  modelplot(list("Constitution" = m_1_twfe,
                           "Law" = m_2_twfe,
                           "Join IO" = m_3_twfe,
                           "Diplomatic mission" = m_4_twfe,
                           "Border patrol" = m_5_twfe,
                           "Justice" = m_6_twfe,
                           "Policing" = m_7_twfe,
                           "Elections" = m_8_twfe,
                           "Currency" = m_9_twfe,
                           "Parallel government" = m_10_twfe,
                           "Local government" = m_11_twfe,
                           "Government in exile" = m_12_twfe,
                           "Tax" = m_13_twfe,
                           "Economic treaty" = m_14_twfe,
                           "Negotiate resource rights" = m_15_twfe,
                           "Education" = m_16_twfe,
                           "Healthcare" = m_17_twfe),
                      coef_rename = c("technologyrebellionIrregular" = "Irregular", 
                                      "technologyrebellionSNC" = "SNC", 
                                      "grptime2" = "Group existance count", 
                                      "rebstrength.num" = "Relative rebel strength",
                                      "e_peinfmor" = "Infant mortality",
                                      "v2exl_legitideol" = "Ideology",
                                      "e_total_resources_income_pc_log" = "Resource income",
                                      "e_pop_log" = "Population",
                                      "e_gdppc_log" = "GDP per capita"),
                      coef_omit = 'Intercept|e_peinfmor|v2exl_legitideol|
          |e_total_resources_income_pc_log|e_gdppc|technologyrebellionSNC|
          |e_pop_log', 
                      facet=T, background = b,
                      draw = FALSE)


# Technology of the Rebellion

f2_data <- modelplot(list("Constitution" = m_a_1_twfe,
                          "Law" = m_a_2_twfe,
                          "Join IO" = m_a_4_twfe,
                          "Diplomatic mission" = m_a_4_twfe,
                          "Border patrol" = m_a_5_twfe,
                          "Justice" = m_a_6_twfe,
                          "Policing" = m_a_7_twfe,
                          "Elections" = m_a_8_twfe,
                          "Currency" = m_a_9_twfe,
                          "Parallel government" = m_a_10_twfe,
                          "Local government" = m_a_11_twfe,
                          "Government in exile" = m_a_12_twfe,
                          "Tax" = m_a_13_twfe,
                          "Economic treaty" = m_a_14_twfe,
                          "Negotiate resource rights" = m_a_15_twfe,
                          "Education" = m_a_16_twfe,
                          "Healthcare" = m_a_17_twfe), 
                     coef_rename = c("technologyrebellionIrregular" = "Irregular", 
                                     "technologyrebellionSNC" = "SNC", 
                                     "grptime2" = "Group existance count", 
                                     "rebstrength.num" = "Relative rebel strength",
                                     "e_peinfmor" = "Infant mortality",
                                     "v2exl_legitideol" = "Ideology",
                                     "e_total_resources_income_pc_log" = "Resource income",
                                     "e_pop_log" = "Population",
                                     "e_gdppc_log" = "GDP per capita"),
                     coef_omit = 'Intercept|e_peinfmor|v2exl_legitideol|
          |e_total_resources_income_pc_log|e_gdppc|technologyrebellionSNC|
          |e_pop_log', 
                     facet=T, background = b, draw = F) 


# plot
f1 <- ggplot(f1_data, aes(
  y = model, x = estimate,
  xmin = conf.low, xmax = conf.high)) +
  ggdist::geom_pointinterval() +
  geom_vline(xintercept = 0, color = 'red', linetype = "longdash",  alpha = 0.5) +
  xlab("Estimate") +
  ylab("Dependent Variable") 

# plot
f2 <- ggplot(f2_data, aes(
  y = model, x = estimate,
  xmin = conf.low, xmax = conf.high)) +
  ggdist::geom_pointinterval() +
  geom_vline(xintercept = 0, color = 'red', linetype = "longdash", alpha = 0.5) +
  xlab("Estimate") +
  ylab("") +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank())

f1 + f2
ggsave("figures/Study_1/Figure_4_ColdWar.png", width = 15, height = 10, units = "cm", dpi = 600)



################################################################################

# Appendix: Regression Models for time after 1991

################################################################################

# Fixed Effects - Year
# Rebel institutions and relative rebel strength

dv <- c("Constitution.count.1", "Law.count.1", "JoinIO.count.1", "DipMiss.count.1", "BorderP.count.1", "Justice.count.1", 
        "Police.count.1", "Elections.count.1", "Currency.count.1", "ParallelGov.count.1", "LocalGov.count.1", 
        "Gov_in_Exile.count.1", "Tax.count.1", "EconTreaty.count.1", "NegResource.count.1", "Edu.count.1", "Health.count.1")

for(i in 1:length(dv)){
  
  model <- paste("m",i, "twfe", sep="_")
  
  m <- feols(as.formula(paste(dv[i], "~ rebstrength.num + e_gdppc + 
                  e_pop_log + e_total_resources_income_pc_log + v2exl_legitideol + 
                              e_peinfmor | Location + year")), cluster = "Location",
             data = subset(qsi.vdem, year %in% c(1992:2020)))
  
  assign(model,m)}

# Rebel institutions and technology of rebellion

for(i in 1:length(dv)){
  
  model <- paste("m_a",i, "twfe", sep="_")
  
  m <- feols(as.formula(paste(dv[i], "~ technologyrebellion + e_gdppc + 
                  e_pop_log + e_total_resources_income_pc_log + v2exl_legitideol + 
                              e_peinfmor | Location + year")), cluster = "Location",
             data = subset(qsi.vdem, year %in% c(1992:2020)))
  
  assign(model,m)}

modelsummary::modelsummary(list(m_1_twfe, m_2_twfe, m_3_twfe, m_4_twfe, m_5_twfe, m_6_twfe, m_7_twfe, m_8_twfe, m_9_twfe, m_10_twfe, 
                                m_11_twfe, m_12_twfe, m_13_twfe, m_14_twfe, m_15_twfe, m_16_twfe, m_17_twfe), 
                           output = "tables/Study_1/Table_A4_after_1991.docx", 
                           estimate = "{estimate}{stars}", 
                           statistic = 'std.error',
                           coef_rename = c("rebstrength.num" = "Rebel Strength", 
                                           "e_gdppc" = "GDP pc (log10)", 
                                           "e_pop_log" = "Pop (log10)", 
                                           "e_total_resources_income_pc_log" = "Ressource Income (log10)", 
                                           "v2exl_legitideol" = "Ideology", 
                                           "e_peinfmor" = "Infant Mortality"))



# Figures

# Rebel strength

f1_data <-  modelplot(list("Constitution" = m_1_twfe,
                           "Law" = m_2_twfe,
                           "Join IO" = m_3_twfe,
                           "Diplomatic mission" = m_4_twfe,
                           "Border patrol" = m_5_twfe,
                           "Justice" = m_6_twfe,
                           "Policing" = m_7_twfe,
                           "Elections" = m_8_twfe,
                           "Currency" = m_9_twfe,
                           "Parallel government" = m_10_twfe,
                           "Local government" = m_11_twfe,
                           "Government in exile" = m_12_twfe,
                           "Tax" = m_13_twfe,
                           "Economic treaty" = m_14_twfe,
                           "Negotiate resource rights" = m_15_twfe,
                           "Education" = m_16_twfe,
                           "Healthcare" = m_17_twfe),
                      coef_rename = c("technologyrebellionIrregular" = "Irregular", 
                                      "technologyrebellionSNC" = "SNC", 
                                      "grptime2" = "Group existance count", 
                                      "rebstrength.num" = "Relative rebel strength",
                                      "e_peinfmor" = "Infant mortality",
                                      "v2exl_legitideol" = "Ideology",
                                      "e_total_resources_income_pc_log" = "Resource income",
                                      "e_pop_log" = "Population",
                                      "e_gdppc_log" = "GDP per capita"),
                      coef_omit = 'Intercept|e_peinfmor|v2exl_legitideol|
          |e_total_resources_income_pc_log|e_gdppc|technologyrebellionSNC|
          |e_pop_log', 
                      facet=T, background = b,
                      draw = FALSE)


# Technology of the Rebellion

f2_data <- modelplot(list("Constitution" = m_a_1_twfe,
                          "Law" = m_a_2_twfe,
                          "Join IO" = m_a_4_twfe,
                          "Diplomatic mission" = m_a_4_twfe,
                          "Border patrol" = m_a_5_twfe,
                          "Justice" = m_a_6_twfe,
                          "Policing" = m_a_7_twfe,
                          "Elections" = m_a_8_twfe,
                          "Currency" = m_a_9_twfe,
                          "Parallel government" = m_a_10_twfe,
                          "Local government" = m_a_11_twfe,
                          "Government in exile" = m_a_12_twfe,
                          "Tax" = m_a_13_twfe,
                          "Economic treaty" = m_a_14_twfe,
                          "Negotiate resource rights" = m_a_15_twfe,
                          "Education" = m_a_16_twfe,
                          "Healthcare" = m_a_17_twfe), 
                     coef_rename = c("technologyrebellionIrregular" = "Irregular", 
                                     "technologyrebellionSNC" = "SNC", 
                                     "grptime2" = "Group existance count", 
                                     "rebstrength.num" = "Relative rebel strength",
                                     "e_peinfmor" = "Infant mortality",
                                     "v2exl_legitideol" = "Ideology",
                                     "e_total_resources_income_pc_log" = "Resource income",
                                     "e_pop_log" = "Population",
                                     "e_gdppc_log" = "GDP per capita"),
                     coef_omit = 'Intercept|e_peinfmor|v2exl_legitideol|
          |e_total_resources_income_pc_log|e_gdppc|technologyrebellionSNC|
          |e_pop_log', 
                     facet=T, background = b, draw = F) 


# plot
f1 <- ggplot(f1_data, aes(
  y = model, x = estimate,
  xmin = conf.low, xmax = conf.high)) +
  ggdist::geom_pointinterval() +
  geom_vline(xintercept = 0, color = 'red', linetype = "longdash",  alpha = 0.5) +
  xlab("Estimate") +
  ylab("Dependent Variable") 

# plot
f2 <- ggplot(f2_data, aes(
  y = model, x = estimate,
  xmin = conf.low, xmax = conf.high)) +
  ggdist::geom_pointinterval() +
  geom_vline(xintercept = 0, color = 'red', linetype = "longdash", alpha = 0.5) +
  xlab("Estimate") +
  ylab("") +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank())

f1 + f2
ggsave("figures/Study_1/Figure_5_after_1991.png", width = 15, height = 10, units = "cm", dpi = 600)



################################################################################

# Appendix: Multiple Imputation for missings 

################################################################################

library(Amelia)



dv <- c("Constitution.count.1", "Law.count.1", "JoinIO.count.1", "DipMiss.count.1", "BorderP.count.1", "Justice.count.1", 
        "Police.count.1", "Elections.count.1", "Currency.count.1", "ParallelGov.count.1", "LocalGov.count.1", 
        "Gov_in_Exile.count.1", "Tax.count.1", "EconTreaty.count.1", "NegResource.count.1", "Edu.count.1", "Health.count.1")

qsi.vdem_amelia <- qsi.vdem %>%
  dplyr::select(dv, rebstrength.num, #technologyrebellion, 
                e_gdppc, e_pop_log, e_total_resources_income_pc_log, v2exl_legitideol, e_peinfmor, Location, year, region)

qsi.vdem_amelia_imputed <- amelia(qsi.vdem_amelia, m = 10, idvars = c(ts = "year", cs= c("Location", "region")))

all_imputations <- bind_rows(unclass(qsi.vdem_amelia_imputed$imputations), .id = "m") %>%
  group_by(m) %>%
  nest()


################################################################################
##### Regression analyses for the rebel strength variables #####


for(i in 1:length(dv)){
  
  model <- paste("m",i, "twfe", sep="_")
  
  m <- all_imputations %>%
    mutate(model = data %>% map(~feols(as.formula(paste(dv[i], "~ rebstrength.num + e_gdppc + 
                  e_pop_log + e_total_resources_income_pc_log + v2exl_legitideol + 
                              e_peinfmor | Location + year")), cluster = "Location",
                                       data = .)), 
           tidied = model %>% map(~ tidy(., conf.int = TRUE)),
           glance = model %>% map(~ glance(.)))    
 
  assign(model,m)
  }

models <- c("m_1_twfe", "m_2_twfe", "m_3_twfe", "m_4_twfe", "m_5_twfe", "m_6_twfe", "m_7_twfe", 
            "m_8_twfe", "m_9_twfe", "m_10_twfe", "m_11_twfe", "m_12_twfe", "m_13_twfe", "m_14_twfe", 
            "m_15_twfe", "m_16_twfe", "m_17_twfe")

params_function <- function(model_name) {
  
  m <- model_name %>%
    unnest(tidied) %>%
    select(m, term, estimate, std.error) %>%
    gather(key, value, estimate, std.error) %>%
    spread(term, value) %>% 
    ungroup()
  return(m)
  
}


m_1_twfe_params <- params_function(m_1_twfe)
m_2_twfe_params <- params_function(m_2_twfe)
m_3_twfe_params <- params_function(m_3_twfe)
m_4_twfe_params <- params_function(m_4_twfe)
m_5_twfe_params <- params_function(m_5_twfe)
m_6_twfe_params <- params_function(m_6_twfe)
m_7_twfe_params <- params_function(m_7_twfe)
m_8_twfe_params <- params_function(m_8_twfe)
m_9_twfe_params <- params_function(m_9_twfe)
m_10_twfe_params <- params_function(m_10_twfe)
m_11_twfe_params <- params_function(m_11_twfe)
m_12_twfe_params <- params_function(m_12_twfe)
m_13_twfe_params <- params_function(m_13_twfe)
m_14_twfe_params <- params_function(m_14_twfe)
m_15_twfe_params <- params_function(m_15_twfe)
m_16_twfe_params <- params_function(m_16_twfe)
m_17_twfe_params <- params_function(m_17_twfe)

coefs_melded_function <- function(model_name) {
  
  model_degree_freedom <- 4566
  
  just_coefs <- model_name %>%
    filter(key == "estimate") %>%
    select(-m, -key)
  
  just_ses <- model_name %>%
    filter(key == "std.error") %>%
    select(-m, -key)
  
  coefs_melded <- mi.meld(just_coefs, just_ses)

  
  melded_summary <- as.data.frame(cbind(t(coefs_melded$q.mi),
                                         t(coefs_melded$se.mi))) %>%
    magrittr::set_colnames(c("estimate", "std.error")) %>%
    mutate(term = rownames(.)) %>%
    select(term, everything()) %>%
    mutate(statistic = estimate / std.error,
           conf.low = estimate + std.error * qt(0.025, model_degree_freedom),
           conf.high = estimate + std.error * qt(0.975, model_degree_freedom),
           p.value = 2 * pt(abs(statistic), model_degree_freedom, lower.tail = FALSE))
  
  
  return(melded_summary)
  
}


melded_summary_m1 <- coefs_melded_function(m_1_twfe_params)
melded_summary_m2 <- coefs_melded_function(m_2_twfe_params)
melded_summary_m3 <- coefs_melded_function(m_3_twfe_params)
melded_summary_m4 <- coefs_melded_function(m_4_twfe_params)
melded_summary_m5 <- coefs_melded_function(m_5_twfe_params)
melded_summary_m6 <- coefs_melded_function(m_6_twfe_params)
melded_summary_m7 <- coefs_melded_function(m_7_twfe_params)
melded_summary_m8 <- coefs_melded_function(m_8_twfe_params)
melded_summary_m9 <- coefs_melded_function(m_9_twfe_params)
melded_summary_m10 <- coefs_melded_function(m_10_twfe_params)
melded_summary_m11 <- coefs_melded_function(m_11_twfe_params)
melded_summary_m12 <- coefs_melded_function(m_12_twfe_params)
melded_summary_m13 <- coefs_melded_function(m_13_twfe_params)
melded_summary_m14 <- coefs_melded_function(m_14_twfe_params)
melded_summary_m15 <- coefs_melded_function(m_15_twfe_params)
melded_summary_m16 <- coefs_melded_function(m_16_twfe_params)
melded_summary_m17 <- coefs_melded_function(m_17_twfe_params)


melded_summary_m1 <- melded_summary_m1 %>%
  mutate(model="Constitution") %>%
  filter(term == "rebstrength.num")
melded_summary_m2 <- melded_summary_m2 %>%
  mutate(model="Law") %>%
  filter(term == "rebstrength.num")
melded_summary_m3 <- melded_summary_m3 %>%
  mutate(model="Join IO") %>%
  filter(term == "rebstrength.num")
melded_summary_m4 <- melded_summary_m4 %>%
  mutate(model="Diplomatic mission") %>%
  filter(term == "rebstrength.num")
melded_summary_m5 <- melded_summary_m5 %>%
  mutate(model="Border patrol") %>%
  filter(term == "rebstrength.num")
melded_summary_m6 <- melded_summary_m6 %>%
  mutate(model="Justice") %>%
  filter(term == "rebstrength.num")
melded_summary_m7 <- melded_summary_m7 %>%
  mutate(model="Policing") %>%
  filter(term == "rebstrength.num")
melded_summary_m8 <- melded_summary_m8 %>%
  mutate(model="Elections") %>%
  filter(term == "rebstrength.num")
melded_summary_m9 <- melded_summary_m9 %>%
  mutate(model="Currency") %>%
  filter(term == "rebstrength.num")
melded_summary_m10 <- melded_summary_m10 %>%
  mutate(model="Parallel government") %>%
  filter(term == "rebstrength.num")
melded_summary_m11 <- melded_summary_m11 %>%
  mutate(model="Local government") %>%
  filter(term == "rebstrength.num")
melded_summary_m12 <- melded_summary_m12 %>%
  mutate(model="Government in exile") %>%
  filter(term == "rebstrength.num")
melded_summary_m13 <- melded_summary_m13 %>%
  mutate(model="Tax") %>%
  filter(term == "rebstrength.num")
melded_summary_m14 <- melded_summary_m14 %>%
  mutate(model="Economic treaty") %>%
  filter(term == "rebstrength.num")
melded_summary_m15 <- melded_summary_m15 %>%
  mutate(model="Negotiate resource rights") %>%
  filter(term == "rebstrength.num")
melded_summary_m16 <- melded_summary_m16 %>%
  mutate(model="Education") %>%
  filter(term == "rebstrength.num")
melded_summary_m17 <- melded_summary_m17 %>%
  mutate(model="Healthcare") %>%
  filter(term == "rebstrength.num")

models_data <- rbind(melded_summary_m1, melded_summary_m2, melded_summary_m3, melded_summary_m4, melded_summary_m5, 
                  melded_summary_m6, melded_summary_m7, melded_summary_m8, melded_summary_m9, melded_summary_m10, 
                  melded_summary_m11, melded_summary_m12, melded_summary_m13, melded_summary_m14, melded_summary_m15, 
                  melded_summary_m16, melded_summary_m17)

models_data <- models_data %>%
  mutate(model = factor(model, levels=c("Constitution", "Law" ,"Join IO","Diplomatic mission" ,"Border patrol","Justice","Policing",
                                        "Elections","Currency","Parallel government","Local government","Government in exile","Tax",
                                        "Economic treaty","Negotiate resource rights","Education","Healthcare")))


# plot
f1 <- ggplot(models_data, aes(
  y = model, x = estimate,
  xmin = conf.low, xmax = conf.high)) +
  ggdist::geom_pointinterval() +
  scale_y_discrete(limits=rev) +
  geom_vline(xintercept = 0, color = 'red', linetype = "longdash",  alpha = 0.5) +
  xlab("Estimate") +
  ylab("Dependent Variable") +
  theme_minimal()




################################################################################
##### Regression analyses for the rebel technology variable #####


dv <- c("Constitution.count.1", "Law.count.1", "JoinIO.count.1", "DipMiss.count.1", "BorderP.count.1", "Justice.count.1", 
        "Police.count.1", "Elections.count.1", "Currency.count.1", "ParallelGov.count.1", "LocalGov.count.1", 
        "Gov_in_Exile.count.1", "Tax.count.1", "EconTreaty.count.1", "NegResource.count.1", "Edu.count.1", "Health.count.1")

qsi.vdem_amelia <- qsi.vdem %>%
  dplyr::select(dv, rebstrength.num, technologyrebellion, 
                e_gdppc, e_pop_log, e_total_resources_income_pc_log, v2exl_legitideol, e_peinfmor, Location, year, region)

qsi.vdem_amelia_imputed <- amelia(qsi.vdem_amelia, m = 10, noms = c("technologyrebellion"),  idvars = c(ts = "year", cs= c("Location", "region")))

all_imputations <- bind_rows(unclass(qsi.vdem_amelia_imputed$imputations), .id = "m") %>%
  group_by(m) %>%
  nest()



for(i in 1:length(dv)){
  
  model <- paste("m_a",i, "twfe", sep="_")
  
  m <- all_imputations %>%
    mutate(model = data %>% map(~feols(as.formula(paste(dv[i], "~ technologyrebellion + e_gdppc + 
                  e_pop_log + e_total_resources_income_pc_log + v2exl_legitideol + 
                              e_peinfmor | Location + year")), cluster = "Location",
                                       data = .)), 
           tidied = model %>% map(~ tidy(., conf.int = TRUE)),
           glance = model %>% map(~ glance(.)))    
  
  assign(model,m)
}

models <- c("m_a_1_twfe", "m_a_2_twfe", "m_a_3_twfe", "m_a_4_twfe", "m_a_5_twfe", "m_a_6_twfe", "m_a_7_twfe", 
            "m_a_8_twfe", "m_a_9_twfe", "m_a_10_twfe", "m_a_11_twfe", "m_a_12_twfe", "m_a_13_twfe", "m_a_14_twfe", 
            "m_a_15_twfe", "m_a_16_twfe", "m_a_17_twfe")

params_function <- function(model_name) {
  
  m <- model_name %>%
    unnest(tidied) %>%
    select(m, term, estimate, std.error) %>%
    gather(key, value, estimate, std.error) %>%
    spread(term, value) %>% 
    ungroup()
  return(m)
  
}


m_1_twfe_params <- params_function(m_a_1_twfe)
m_2_twfe_params <- params_function(m_a_2_twfe)
m_3_twfe_params <- params_function(m_a_3_twfe)
m_4_twfe_params <- params_function(m_a_4_twfe)
m_5_twfe_params <- params_function(m_a_5_twfe)
m_6_twfe_params <- params_function(m_a_6_twfe)
m_7_twfe_params <- params_function(m_a_7_twfe)
m_8_twfe_params <- params_function(m_a_8_twfe)
m_9_twfe_params <- params_function(m_a_9_twfe)
m_10_twfe_params <- params_function(m_a_10_twfe)
m_11_twfe_params <- params_function(m_a_11_twfe)
m_12_twfe_params <- params_function(m_a_12_twfe)
m_13_twfe_params <- params_function(m_a_13_twfe)
m_14_twfe_params <- params_function(m_a_14_twfe)
m_15_twfe_params <- params_function(m_a_15_twfe)
m_16_twfe_params <- params_function(m_a_16_twfe)
m_17_twfe_params <- params_function(m_a_17_twfe)

coefs_melded_function <- function(model_name) {
  
  model_degree_freedom <- 4566
  
  just_coefs <- model_name %>%
    filter(key == "estimate") %>%
    select(-m, -key)
  
  just_ses <- model_name %>%
    filter(key == "std.error") %>%
    select(-m, -key)
  
  coefs_melded <- mi.meld(just_coefs, just_ses)
  
  
  melded_summary <- as.data.frame(cbind(t(coefs_melded$q.mi),
                                        t(coefs_melded$se.mi))) %>%
    magrittr::set_colnames(c("estimate", "std.error")) %>%
    mutate(term = rownames(.)) %>%
    select(term, everything()) %>%
    mutate(statistic = estimate / std.error,
           conf.low = estimate + std.error * qt(0.025, model_degree_freedom),
           conf.high = estimate + std.error * qt(0.975, model_degree_freedom),
           p.value = 2 * pt(abs(statistic), model_degree_freedom, lower.tail = FALSE))
  
  
  return(melded_summary)
  
}


melded_summary_m1 <- coefs_melded_function(m_1_twfe_params)
melded_summary_m2 <- coefs_melded_function(m_2_twfe_params)
melded_summary_m3 <- coefs_melded_function(m_3_twfe_params)
melded_summary_m4 <- coefs_melded_function(m_4_twfe_params)
melded_summary_m5 <- coefs_melded_function(m_5_twfe_params)
melded_summary_m6 <- coefs_melded_function(m_6_twfe_params)
melded_summary_m7 <- coefs_melded_function(m_7_twfe_params)
melded_summary_m8 <- coefs_melded_function(m_8_twfe_params)
melded_summary_m9 <- coefs_melded_function(m_9_twfe_params)
melded_summary_m10 <- coefs_melded_function(m_10_twfe_params)
melded_summary_m11 <- coefs_melded_function(m_11_twfe_params)
melded_summary_m12 <- coefs_melded_function(m_12_twfe_params)
melded_summary_m13 <- coefs_melded_function(m_13_twfe_params)
melded_summary_m14 <- coefs_melded_function(m_14_twfe_params)
melded_summary_m15 <- coefs_melded_function(m_15_twfe_params)
melded_summary_m16 <- coefs_melded_function(m_16_twfe_params)
melded_summary_m17 <- coefs_melded_function(m_17_twfe_params)


melded_summary_m1 <- melded_summary_m1 %>%
  mutate(model="Constitution") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m2 <- melded_summary_m2 %>%
  mutate(model="Law") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m3 <- melded_summary_m3 %>%
  mutate(model="Join IO") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m4 <- melded_summary_m4 %>%
  mutate(model="Diplomatic mission") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m5 <- melded_summary_m5 %>%
  mutate(model="Border patrol") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m6 <- melded_summary_m6 %>%
  mutate(model="Justice") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m7 <- melded_summary_m7 %>%
  mutate(model="Policing") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m8 <- melded_summary_m8 %>%
  mutate(model="Elections") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m9 <- melded_summary_m9 %>%
  mutate(model="Currency") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m10 <- melded_summary_m10 %>%
  mutate(model="Parallel government") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m11 <- melded_summary_m11 %>%
  mutate(model="Local government") %>%
  filter(term == "rebstrength.num")
melded_summary_m12 <- melded_summary_m12 %>%
  mutate(model="Government in exile") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m13 <- melded_summary_m13 %>%
  mutate(model="Tax") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m14 <- melded_summary_m14 %>%
  mutate(model="Economic treaty") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m15 <- melded_summary_m15 %>%
  mutate(model="Negotiate resource rights") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m16 <- melded_summary_m16 %>%
  mutate(model="Education") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m17 <- melded_summary_m17 %>%
  mutate(model="Healthcare") %>%
  filter(term == "technologyrebellionIrregular")

models_data <- rbind(melded_summary_m1, melded_summary_m2, melded_summary_m3, melded_summary_m4, melded_summary_m5, 
                     melded_summary_m6, melded_summary_m7, melded_summary_m8, melded_summary_m9, melded_summary_m10, 
                     melded_summary_m11, melded_summary_m12, melded_summary_m13, melded_summary_m14, melded_summary_m15, 
                     melded_summary_m16, melded_summary_m17)

models_data <- models_data %>%
  mutate(model = factor(model, levels=c("Constitution", "Law" ,"Join IO","Diplomatic mission" ,"Border patrol","Justice","Policing",
                                        "Elections","Currency","Parallel government","Local government","Government in exile","Tax",
                                        "Economic treaty","Negotiate resource rights","Education","Healthcare")))


# plot
f2 <- ggplot(models_data, aes(
  y = model, x = estimate,
  xmin = conf.low, xmax = conf.high)) +
  ggdist::geom_pointinterval() +
  scale_y_discrete(limits=rev) +
  geom_vline(xintercept = 0, color = 'red', linetype = "longdash",  alpha = 0.5) +
  xlab("Estimate") +
  ylab("") +
  theme_minimal() +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank())


f1 + f2
ggsave("figures/Study_1/Figure_6_MI.png", width = 15, height = 10, units = "cm", dpi = 600)


################################################################################

# Appendix: Multiple Imputation with region FE

################################################################################

dv <- c("Constitution.count.1", "Law.count.1", "JoinIO.count.1", "DipMiss.count.1", "BorderP.count.1", "Justice.count.1", 
        "Police.count.1", "Elections.count.1", "Currency.count.1", "ParallelGov.count.1", "LocalGov.count.1", 
        "Gov_in_Exile.count.1", "Tax.count.1", "EconTreaty.count.1", "NegResource.count.1", "Edu.count.1", "Health.count.1")

qsi.vdem_amelia <- qsi.vdem %>%
  dplyr::select(dv, rebstrength.num, #technologyrebellion, 
                e_gdppc, e_pop_log, e_total_resources_income_pc_log, v2exl_legitideol, e_peinfmor, Location, year, region)

qsi.vdem_amelia_imputed <- amelia(qsi.vdem_amelia, m = 10, idvars = c(ts = "year", cs= c("Location", "region")))

all_imputations <- bind_rows(unclass(qsi.vdem_amelia_imputed$imputations), .id = "m") %>%
  group_by(m) %>%
  nest()


################################################################################
##### Regression analyses for the rebel strength variables #####


for(i in 1:length(dv)){
  
  model <- paste("m",i, "twfe", sep="_")
  
  m <- all_imputations %>%
    mutate(model = data %>% map(~feols(as.formula(paste(dv[i], "~ rebstrength.num + e_gdppc + 
                  e_pop_log + e_total_resources_income_pc_log + v2exl_legitideol + 
                              e_peinfmor | region + year")), cluster = "Location",
                                       data = .)), 
           tidied = model %>% map(~ tidy(., conf.int = TRUE)),
           glance = model %>% map(~ glance(.)))    
  
  assign(model,m)
}

models <- c("m_1_twfe", "m_2_twfe", "m_3_twfe", "m_4_twfe", "m_5_twfe", "m_6_twfe", "m_7_twfe", 
            "m_8_twfe", "m_9_twfe", "m_10_twfe", "m_11_twfe", "m_12_twfe", "m_13_twfe", "m_14_twfe", 
            "m_15_twfe", "m_16_twfe", "m_17_twfe")

params_function <- function(model_name) {
  
  m <- model_name %>%
    unnest(tidied) %>%
    select(m, term, estimate, std.error) %>%
    gather(key, value, estimate, std.error) %>%
    spread(term, value) %>% 
    ungroup()
  return(m)
  
}


m_1_twfe_params <- params_function(m_1_twfe)
m_2_twfe_params <- params_function(m_2_twfe)
m_3_twfe_params <- params_function(m_3_twfe)
m_4_twfe_params <- params_function(m_4_twfe)
m_5_twfe_params <- params_function(m_5_twfe)
m_6_twfe_params <- params_function(m_6_twfe)
m_7_twfe_params <- params_function(m_7_twfe)
m_8_twfe_params <- params_function(m_8_twfe)
m_9_twfe_params <- params_function(m_9_twfe)
m_10_twfe_params <- params_function(m_10_twfe)
m_11_twfe_params <- params_function(m_11_twfe)
m_12_twfe_params <- params_function(m_12_twfe)
m_13_twfe_params <- params_function(m_13_twfe)
m_14_twfe_params <- params_function(m_14_twfe)
m_15_twfe_params <- params_function(m_15_twfe)
m_16_twfe_params <- params_function(m_16_twfe)
m_17_twfe_params <- params_function(m_17_twfe)

coefs_melded_function <- function(model_name) {
  
  model_degree_freedom <- 4566
  
  just_coefs <- model_name %>%
    filter(key == "estimate") %>%
    select(-m, -key)
  
  just_ses <- model_name %>%
    filter(key == "std.error") %>%
    select(-m, -key)
  
  coefs_melded <- mi.meld(just_coefs, just_ses)
  
  
  melded_summary <- as.data.frame(cbind(t(coefs_melded$q.mi),
                                        t(coefs_melded$se.mi))) %>%
    magrittr::set_colnames(c("estimate", "std.error")) %>%
    mutate(term = rownames(.)) %>%
    select(term, everything()) %>%
    mutate(statistic = estimate / std.error,
           conf.low = estimate + std.error * qt(0.025, model_degree_freedom),
           conf.high = estimate + std.error * qt(0.975, model_degree_freedom),
           p.value = 2 * pt(abs(statistic), model_degree_freedom, lower.tail = FALSE))
  
  
  return(melded_summary)
  
}


melded_summary_m1 <- coefs_melded_function(m_1_twfe_params)
melded_summary_m2 <- coefs_melded_function(m_2_twfe_params)
melded_summary_m3 <- coefs_melded_function(m_3_twfe_params)
melded_summary_m4 <- coefs_melded_function(m_4_twfe_params)
melded_summary_m5 <- coefs_melded_function(m_5_twfe_params)
melded_summary_m6 <- coefs_melded_function(m_6_twfe_params)
melded_summary_m7 <- coefs_melded_function(m_7_twfe_params)
melded_summary_m8 <- coefs_melded_function(m_8_twfe_params)
melded_summary_m9 <- coefs_melded_function(m_9_twfe_params)
melded_summary_m10 <- coefs_melded_function(m_10_twfe_params)
melded_summary_m11 <- coefs_melded_function(m_11_twfe_params)
melded_summary_m12 <- coefs_melded_function(m_12_twfe_params)
melded_summary_m13 <- coefs_melded_function(m_13_twfe_params)
melded_summary_m14 <- coefs_melded_function(m_14_twfe_params)
melded_summary_m15 <- coefs_melded_function(m_15_twfe_params)
melded_summary_m16 <- coefs_melded_function(m_16_twfe_params)
melded_summary_m17 <- coefs_melded_function(m_17_twfe_params)


melded_summary_m1 <- melded_summary_m1 %>%
  mutate(model="Constitution") %>%
  filter(term == "rebstrength.num")
melded_summary_m2 <- melded_summary_m2 %>%
  mutate(model="Law") %>%
  filter(term == "rebstrength.num")
melded_summary_m3 <- melded_summary_m3 %>%
  mutate(model="Join IO") %>%
  filter(term == "rebstrength.num")
melded_summary_m4 <- melded_summary_m4 %>%
  mutate(model="Diplomatic mission") %>%
  filter(term == "rebstrength.num")
melded_summary_m5 <- melded_summary_m5 %>%
  mutate(model="Border patrol") %>%
  filter(term == "rebstrength.num")
melded_summary_m6 <- melded_summary_m6 %>%
  mutate(model="Justice") %>%
  filter(term == "rebstrength.num")
melded_summary_m7 <- melded_summary_m7 %>%
  mutate(model="Policing") %>%
  filter(term == "rebstrength.num")
melded_summary_m8 <- melded_summary_m8 %>%
  mutate(model="Elections") %>%
  filter(term == "rebstrength.num")
melded_summary_m9 <- melded_summary_m9 %>%
  mutate(model="Currency") %>%
  filter(term == "rebstrength.num")
melded_summary_m10 <- melded_summary_m10 %>%
  mutate(model="Parallel government") %>%
  filter(term == "rebstrength.num")
melded_summary_m11 <- melded_summary_m11 %>%
  mutate(model="Local government") %>%
  filter(term == "rebstrength.num")
melded_summary_m12 <- melded_summary_m12 %>%
  mutate(model="Government in exile") %>%
  filter(term == "rebstrength.num")
melded_summary_m13 <- melded_summary_m13 %>%
  mutate(model="Tax") %>%
  filter(term == "rebstrength.num")
melded_summary_m14 <- melded_summary_m14 %>%
  mutate(model="Economic treaty") %>%
  filter(term == "rebstrength.num")
melded_summary_m15 <- melded_summary_m15 %>%
  mutate(model="Negotiate resource rights") %>%
  filter(term == "rebstrength.num")
melded_summary_m16 <- melded_summary_m16 %>%
  mutate(model="Education") %>%
  filter(term == "rebstrength.num")
melded_summary_m17 <- melded_summary_m17 %>%
  mutate(model="Healthcare") %>%
  filter(term == "rebstrength.num")

models_data <- rbind(melded_summary_m1, melded_summary_m2, melded_summary_m3, melded_summary_m4, melded_summary_m5, 
                     melded_summary_m6, melded_summary_m7, melded_summary_m8, melded_summary_m9, melded_summary_m10, 
                     melded_summary_m11, melded_summary_m12, melded_summary_m13, melded_summary_m14, melded_summary_m15, 
                     melded_summary_m16, melded_summary_m17)

models_data <- models_data %>%
  mutate(model = factor(model, levels=c("Constitution", "Law" ,"Join IO","Diplomatic mission" ,"Border patrol","Justice","Policing",
                                        "Elections","Currency","Parallel government","Local government","Government in exile","Tax",
                                        "Economic treaty","Negotiate resource rights","Education","Healthcare")))


# plot
f1 <- ggplot(models_data, aes(
  y = model, x = estimate,
  xmin = conf.low, xmax = conf.high)) +
  ggdist::geom_pointinterval() +
  scale_y_discrete(limits=rev) +
  geom_vline(xintercept = 0, color = 'red', linetype = "longdash",  alpha = 0.5) +
  xlab("Estimate") +
  ylab("Dependent Variable") +
  theme_minimal()




################################################################################
##### Regression analyses for the rebel technology variable #####


dv <- c("Constitution.count.1", "Law.count.1", "JoinIO.count.1", "DipMiss.count.1", "BorderP.count.1", "Justice.count.1", 
        "Police.count.1", "Elections.count.1", "Currency.count.1", "ParallelGov.count.1", "LocalGov.count.1", 
        "Gov_in_Exile.count.1", "Tax.count.1", "EconTreaty.count.1", "NegResource.count.1", "Edu.count.1", "Health.count.1")

qsi.vdem_amelia <- qsi.vdem %>%
  dplyr::select(dv, rebstrength.num, technologyrebellion, 
                e_gdppc, e_pop_log, e_total_resources_income_pc_log, v2exl_legitideol, e_peinfmor, Location, year, region)

qsi.vdem_amelia_imputed <- amelia(qsi.vdem_amelia, m = 10, noms = c("technologyrebellion"),  idvars = c(ts = "year", cs= c("Location", "region")))

all_imputations <- bind_rows(unclass(qsi.vdem_amelia_imputed$imputations), .id = "m") %>%
  group_by(m) %>%
  nest()



for(i in 1:length(dv)){
  
  model <- paste("m_a",i, "twfe", sep="_")
  
  m <- all_imputations %>%
    mutate(model = data %>% map(~feols(as.formula(paste(dv[i], "~ technologyrebellion + e_gdppc + 
                  e_pop_log + e_total_resources_income_pc_log + v2exl_legitideol + 
                              e_peinfmor | region + year")), cluster = "Location",
                                       data = .)), 
           tidied = model %>% map(~ tidy(., conf.int = TRUE)),
           glance = model %>% map(~ glance(.)))    
  
  assign(model,m)
}

models <- c("m_a_1_twfe", "m_a_2_twfe", "m_a_3_twfe", "m_a_4_twfe", "m_a_5_twfe", "m_a_6_twfe", "m_a_7_twfe", 
            "m_a_8_twfe", "m_a_9_twfe", "m_a_10_twfe", "m_a_11_twfe", "m_a_12_twfe", "m_a_13_twfe", "m_a_14_twfe", 
            "m_a_15_twfe", "m_a_16_twfe", "m_a_17_twfe")

params_function <- function(model_name) {
  
  m <- model_name %>%
    unnest(tidied) %>%
    select(m, term, estimate, std.error) %>%
    gather(key, value, estimate, std.error) %>%
    spread(term, value) %>% 
    ungroup()
  return(m)
  
}


m_1_twfe_params <- params_function(m_a_1_twfe)
m_2_twfe_params <- params_function(m_a_2_twfe)
m_3_twfe_params <- params_function(m_a_3_twfe)
m_4_twfe_params <- params_function(m_a_4_twfe)
m_5_twfe_params <- params_function(m_a_5_twfe)
m_6_twfe_params <- params_function(m_a_6_twfe)
m_7_twfe_params <- params_function(m_a_7_twfe)
m_8_twfe_params <- params_function(m_a_8_twfe)
m_9_twfe_params <- params_function(m_a_9_twfe)
m_10_twfe_params <- params_function(m_a_10_twfe)
m_11_twfe_params <- params_function(m_a_11_twfe)
m_12_twfe_params <- params_function(m_a_12_twfe)
m_13_twfe_params <- params_function(m_a_13_twfe)
m_14_twfe_params <- params_function(m_a_14_twfe)
m_15_twfe_params <- params_function(m_a_15_twfe)
m_16_twfe_params <- params_function(m_a_16_twfe)
m_17_twfe_params <- params_function(m_a_17_twfe)

coefs_melded_function <- function(model_name) {
  
  model_degree_freedom <- 4566
  
  just_coefs <- model_name %>%
    filter(key == "estimate") %>%
    select(-m, -key)
  
  just_ses <- model_name %>%
    filter(key == "std.error") %>%
    select(-m, -key)
  
  coefs_melded <- mi.meld(just_coefs, just_ses)
  
  
  melded_summary <- as.data.frame(cbind(t(coefs_melded$q.mi),
                                        t(coefs_melded$se.mi))) %>%
    magrittr::set_colnames(c("estimate", "std.error")) %>%
    mutate(term = rownames(.)) %>%
    select(term, everything()) %>%
    mutate(statistic = estimate / std.error,
           conf.low = estimate + std.error * qt(0.025, model_degree_freedom),
           conf.high = estimate + std.error * qt(0.975, model_degree_freedom),
           p.value = 2 * pt(abs(statistic), model_degree_freedom, lower.tail = FALSE))
  
  
  return(melded_summary)
  
}


melded_summary_m1 <- coefs_melded_function(m_1_twfe_params)
melded_summary_m2 <- coefs_melded_function(m_2_twfe_params)
melded_summary_m3 <- coefs_melded_function(m_3_twfe_params)
melded_summary_m4 <- coefs_melded_function(m_4_twfe_params)
melded_summary_m5 <- coefs_melded_function(m_5_twfe_params)
melded_summary_m6 <- coefs_melded_function(m_6_twfe_params)
melded_summary_m7 <- coefs_melded_function(m_7_twfe_params)
melded_summary_m8 <- coefs_melded_function(m_8_twfe_params)
melded_summary_m9 <- coefs_melded_function(m_9_twfe_params)
melded_summary_m10 <- coefs_melded_function(m_10_twfe_params)
melded_summary_m11 <- coefs_melded_function(m_11_twfe_params)
melded_summary_m12 <- coefs_melded_function(m_12_twfe_params)
melded_summary_m13 <- coefs_melded_function(m_13_twfe_params)
melded_summary_m14 <- coefs_melded_function(m_14_twfe_params)
melded_summary_m15 <- coefs_melded_function(m_15_twfe_params)
melded_summary_m16 <- coefs_melded_function(m_16_twfe_params)
melded_summary_m17 <- coefs_melded_function(m_17_twfe_params)


melded_summary_m1 <- melded_summary_m1 %>%
  mutate(model="Constitution") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m2 <- melded_summary_m2 %>%
  mutate(model="Law") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m3 <- melded_summary_m3 %>%
  mutate(model="Join IO") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m4 <- melded_summary_m4 %>%
  mutate(model="Diplomatic mission") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m5 <- melded_summary_m5 %>%
  mutate(model="Border patrol") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m6 <- melded_summary_m6 %>%
  mutate(model="Justice") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m7 <- melded_summary_m7 %>%
  mutate(model="Policing") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m8 <- melded_summary_m8 %>%
  mutate(model="Elections") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m9 <- melded_summary_m9 %>%
  mutate(model="Currency") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m10 <- melded_summary_m10 %>%
  mutate(model="Parallel government") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m11 <- melded_summary_m11 %>%
  mutate(model="Local government") %>%
  filter(term == "rebstrength.num")
melded_summary_m12 <- melded_summary_m12 %>%
  mutate(model="Government in exile") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m13 <- melded_summary_m13 %>%
  mutate(model="Tax") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m14 <- melded_summary_m14 %>%
  mutate(model="Economic treaty") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m15 <- melded_summary_m15 %>%
  mutate(model="Negotiate resource rights") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m16 <- melded_summary_m16 %>%
  mutate(model="Education") %>%
  filter(term == "technologyrebellionIrregular")
melded_summary_m17 <- melded_summary_m17 %>%
  mutate(model="Healthcare") %>%
  filter(term == "technologyrebellionIrregular")

models_data <- rbind(melded_summary_m1, melded_summary_m2, melded_summary_m3, melded_summary_m4, melded_summary_m5, 
                     melded_summary_m6, melded_summary_m7, melded_summary_m8, melded_summary_m9, melded_summary_m10, 
                     melded_summary_m11, melded_summary_m12, melded_summary_m13, melded_summary_m14, melded_summary_m15, 
                     melded_summary_m16, melded_summary_m17)

models_data <- models_data %>%
  mutate(model = factor(model, levels=c("Constitution", "Law" ,"Join IO","Diplomatic mission" ,"Border patrol","Justice","Policing",
                                        "Elections","Currency","Parallel government","Local government","Government in exile","Tax",
                                        "Economic treaty","Negotiate resource rights","Education","Healthcare")))


# plot
f2 <- ggplot(models_data, aes(
  y = model, x = estimate,
  xmin = conf.low, xmax = conf.high)) +
  ggdist::geom_pointinterval() +
  scale_y_discrete(limits=rev) +
  geom_vline(xintercept = 0, color = 'red', linetype = "longdash",  alpha = 0.5) +
  xlab("Estimate") +
  ylab("") +
  theme_minimal() +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank())


f1 + f2
ggsave("figures/Study_1/Figure_7_MI_Region_FE.png", width = 15, height = 10, units = "cm", dpi = 600)

